Reading data


Textinho explicando o que foi feito no primeiro notebook

[1]:
FILE_PATH = '/content/drive/MyDrive/01 - Iniciação Científica/IC-CoRoT_Kepler/resampled_files/RESAMPLED_EN2_STAR_CHR_0101086161_20070516T060226_20071005T074409.csv'
[2]:
import pandas as pd

data = pd.read_csv(FILE_PATH)
data.head()
[2]:
DATE WHITEFLUX
0 54236.757582 112521.329834
1 54236.767013 112758.045853
2 54236.776445 112943.042225
3 54236.785876 112562.266242
4 54236.795308 112789.303079
[3]:
import numpy as np

time = data.DATE.to_numpy()
flux = data.WHITEFLUX.to_numpy()
[4]:
from utils import *
[5]:
curve = lightcurve.LightCurve(time, flux)
curve.plot()

Ideal Lowpass Filter - Traduzir


Um filtro bi-dimensional passa-baixa que deixa passar todas as frequências em um círculo de raio \(D_0\) a partir da origem e remove todas as frequências fora desse círculo é chamado de filtro passa-baixa ideal (ILPF) e é descrito como

\[\begin{split} H(u) = \begin{cases} 1, &\text{if } D(u) \le D_0 \\ 0, &\text{if } D(u) \ge D_0 \end{cases}\end{split}\]

onde \(D_0\) é uma constante positiva, e \(D(u)\) é a distância entre um ponto \(u\) até o centro do retângulo de frequência, ou seja, é definido por

\[D(u) = (u-P/2) \qquad (1)\]

sendo \(P\) o tamanho do vetor original preenchido (padded).

O ponto de transição entre \(H(u) = 1\) e \(H(u) = 0\) é chamado de frequência de corte

[ ]:
def ideal_filter(array, fourier_transform, cutoff_freq):
  n_time = len(array)
  D0 = cutoff_freq * n_time
  xc = n_time

  for i in range(len(fourier_transform)):
    if fourier_transform[i] > D0:
      fourier_transform[i] = 0

  y_filtered = np.real(np.fft.ifft(fourier_transform))
  y_filtered += (array.mean() - y_filtered.mean())

  return y_filtered
[ ]:
def ideal_filter_result(array, cutoff_freq, numExpansion):
  Filter = filters.FrequencyDomainFiltering()

  Filter.expand_borders(y, numExpansion)
  y_expanded = Filter.getExpandedBorders

  y_filtered = ideal_filter(y_expanded, np.fft.fft(y_expanded), cutoff_freq)

  Filter.remove_expanded_borders(y_filtered, numExpansion)
  y_filtered = Filter.getNoExpanded

  return y_filtered

Plotting some results

[ ]:
filtered = curve.ideal_lowpass_filter(0.1)
filtered.view_filter_results()
[ ]:
filtered = curve.ideal_lowpass_filter(0.6)
filtered.view_filter_results()

Saving filtered data

[ ]:
cutoff_freqs = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]
[ ]:
PATH_DIR = 'C:/Users/guisa/Google Drive/01 - Iniciação Científica/02 - Datasets/exoplanets_confirmed/filters/ideal/cutoffFreq0'
[ ]:
%time

for cutoff_freq in cutoff_freqs:
   # Saving filtered data
   # filters.export_results_csv(PATH_DIR=PATH_DIR+str(int(cutoff_freq*10)), filter_technique='ideal', cutoff_freq=cutoff_freq, order=None, numNei=None)
CPU times: user 3 µs, sys: 0 ns, total: 3 µs
Wall time: 7.15 µs
[ ]:

[ ]:

[ ]:

Gaussian Lowpass Filter


The transfer function of a Gaussian 1-D lowpass filter (GLPFs) is defined by

\[H(u) = e^{-D^2(u)/2D^2_0}\]

where \(D(u)\) and \(D_0\) was defined on Eq. (1).

Note. The cutoff frequency must be given in Nyquist.

[ ]:
def gaussian_array(array, fourier_transform, cutoff_freq):
  # Extrating information of the signal
  n_time = len(array)
  D0 = cutoff_freq * n_time
  xc = n_time

  # Creating the filter array
  len_filter = len(fourier_transform)
  filter = np.zeros(len_filter)

  for i in range(len_filter):
    filter[i] = exp( (-(i-(xc-1.0))**2)/(2*((D0 * n_time)**2)) )

  return filter

Plotting some results - Completar

Note. It was observed that, at low cutoff frequencies (0.1 and 0.2 Nyquist), there is an effect in which the filtered curve undergoes a vertical displacement, as if a constant were added to the curve. This will not affect the modeling of the curves, but to prevent this effect from causing disturbances in the visualization of the results, the following operation is perfomed:

\[Flux\space Filtered\space += [mean(Raw\space Flux) - mean(Flux\space Filtered)]\]
[ ]:
filtered = curve.gaussian_lowpass_filter(cutoff_freq=0.1)
filtered.view_filter_results()
[ ]:
filtered = curve.gaussian_lowpass_filter(cutoff_freq=0.4)
filtered.view_filter_results()

Saving filtered data

[ ]:
cutoff_freqs = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]
[ ]:
PATH_DIR = 'C:/Users/guisa/Google Drive/01 - Iniciação Científica/02 - Datasets/exoplanets_confirmed/filters/gaussian/cutoffFreq0'
[ ]:
%time

for cutoff_freq in cutoff_freqs:
   # Saving filtered data
   # filters.export_results_csv(PATH_DIR=PATH_DIR+str(int(cutoff_freq*10)), filter_technique='gaussian', cutoff_freq=cutoff_freq, order=None, numNei=None)
CPU times: user 3 µs, sys: 0 ns, total: 3 µs
Wall time: 8.34 µs

Butterworth Lowpass Filter


The transfer function of a Butterworth 1-D lowpass filter (BLPF) of order \(n\), and with cutoff frequency at a distance \(D_0\) from the origin, is defined as

\[H(u) = \frac{1}{ 1 + [D(u) / D_0]^{2n} }\]

where \(D(u)\) and \(D_0\) was defined on Eq. (1).

By the definition, the Butterworth filter have two free parameters: the cutoff frequency and the filtering order. Then, we can modify both, as we can see on the code cell below, intending to have the best results possibles.

Note. The cutoff frequency must be given in Nyquist.

[ ]:
def butterworth_array(array, fourier_transform, cutoff_freq, order):
  # Extrating information of the signal
  n_time = len(array)
  D0 = cutoff_freq * n_time
  xc = n_time

  # Creating the filter array
  len_filter = len(fourier_transform)
  filter = np.zeros(len_filter)

  for i in range(len_filter):
    filter[i] = 1.0 / (1.0+(abs(i-(xc-1.0))/D0)**(2.0*order))

  return filter

Plotting some results

[ ]:
filtered = curve.butterworth_lowpass_filter(2, 0.1)
filtered.view_filter_results()
[ ]:
filtered = curve.butterworth_lowpass_filter(4, 0.3)
filtered.view_filter_results()

Saving filtered data

[ ]:
orders = [1, 2, 3, 4, 5, 6]
cutoff_freqs = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]
[ ]:
%time

for order in orders:
  for cutoff_freq in cutoff_freqs:
    # Saving filtered data
    PATH_DIR = 'C:/Users/guisa/Google Drive/01 - Iniciação Científica/02 - Datasets/exoplanets_confirmed/filters/butterworth/order'+str(order)+'/cutoffFreq0'+str(int(cutoff_freq*10))
    # filters.export_results_csv(PATH_DIR=PATH_DIR+str(int(cutoff_freq*10)), filter_technique='butterworth', cutoff_freq=cutoff_freq, order=orders, numNei=None)
CPU times: user 3 µs, sys: 0 ns, total: 3 µs
Wall time: 6.44 µs

Bessel Lowpass Filter


The transfer function, \(H(s)\), of a Bessel lowpass filter is defined by

\[\begin{split}H(s) = \frac{\theta_n (0) }{\theta_n (s/ \omega_0)} \\\end{split}\]

where

\[\theta_n (s) = \sum_{k=0}^{n} a_k s^k\]

and

\[\begin{split}a_k = \frac{(2n-k)!}{2^{n-k}k!(n-k)!} \\\end{split}\]

Showing how it works…

Parameters

[ ]:
order = 2
cutoff_freq = 0.6

Control lib

[ ]:
from control import *
[ ]:
### Computing ak

from math import factorial

coef = []
i = 0
while i <= order:
  ak = (factorial(2*order - i)) / ( 2**(order - i)*factorial(i)*factorial(order - i) )
  # print(ak)
  coef.append(ak)
  i += 1

print(coef)
[3.0, 3.0, 1.0]
[ ]:
### Computing θn(s)

s = TransferFunction.s
theta_array = []
k = 0
for k in range(order+1):
  theta_n = coef[k] * (s**k)
  theta_array.append(theta_n)
  # numerical_numerator = coef[0]
  # print(theta_n)

print(theta_array[0])
print(theta_array[1])
print(theta_array[2])

3
-
1


3 s
---
 1


s^2
---
 1

[ ]:
### Computing H(s)

coef_numerator = theta_array[0]

list_denominator = theta_array[:]
[ ]:
denominator = 0
for item in list_denominator:
  denominator += item

print(denominator)

s^2 + 3 s + 3
-------------
      1

[ ]:
### Filling in transfer function

G = coef_numerator / denominator

print(G)
print(type(G))

      3
-------------
s^2 + 3 s + 3

<class 'control.xferfcn.TransferFunction'>

Applying

[ ]:
def bessel(array, fourier_transform, cutoff_freq, order):

  # Extracting features from signal
  n_time = len(array)
  D0 = cutoff_freq * n_time
  xc = n_time

  # Creating the bessel transfer function array
  len_filter = len(fourier_transform)
  filter = np.zeros(len_filter)
  i=0

  for i in range(len_filter):
    filter[i] = np.real(evalfr(G, ( np.abs(i-(xc-1.0))/D0 )))

  return filter

Plotting some results

[ ]:
filtered = curve.bessel_lowpass_filter(2, 0.1, numExpansion=100)
filtered.view_filter_results()
[ ]:
filtered = curve.bessel_lowpass_filter(4, 0.3, numExpansion=100)
filtered.view_filter_results()

Saving filtered data

[ ]:
orders = [1, 2, 3, 4, 5, 6]
cutoff_freqs = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]
[ ]:
%time

for order in orders:
  for cutoff_freq in cutoff_freqs:
    # Saving filtered data
    PATH_DIR = 'C:/Users/guisa/Google Drive/01 - Iniciação Científica/02 - Datasets/exoplanets_confirmed/filters/bessel/order'+str(order)+'/cutoffFreq0'+str(int(cutoff_freq*10))
    # filters.export_results_csv(PATH_DIR=PATH_DIR+str(int(cutoff_freq*10)), filter_technique='butterworth', cutoff_freq=cutoff_freq, order=orders, numNei=None)
CPU times: user 3 µs, sys: 1 µs, total: 4 µs
Wall time: 5.72 µs

Median Filter


O filtro mediano, por sua vez, é aplicado de uma forma consideravelmente diferente, na qual cada valor dos dados filtrados corresponde a uma mediana de um grupo de valores adjacentes nos dados originais. Esse filtro já é utilizado com certa frequência em análises de curvas de luz.

[ ]:
from scipy.signal import medfilt

def median_filter(array, window_size):
    return medfilt(array, window_size)

Reference: Scipy Documentation

Plotting some results

[ ]:
filtered = curve.median_filter(3)
filtered.view_filter_results()
[ ]:
filtered = curve.median_filter(9)
filtered.view_filter_results()

Saving filtered data

[ ]:
neighbors = [3, 5, 7, 9, 11]
[ ]:
%time

for neighbor in neighbors:
  # Saving filtered data
  PATH_DIR = 'C:/Users/guisa/Google Drive/01 - Iniciação Científica/02 - Datasets/exoplanets_confirmed/filters/median/numNei'+str(int(neighbor))
  #filters.export_results_csv(PATH_DIR, 'median', cutoff_freq=None, order=None, numNei=neighbor)
CPU times: user 4 µs, sys: 1e+03 ns, total: 5 µs
Wall time: 8.34 µs